{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib as mpl\n",
    "import matplotlib.colors\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from sklearn.mixture import GaussianMixture\n",
    "from sklearn.metrics.pairwise import pairwise_distances_argmin"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# 解决中文显示问题\n",
    "mpl.rcParams['font.sans-serif'] = [u'SimHei']\n",
    "mpl.rcParams['axes.unicode_minus'] = False"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# 设置在jupyter中matplotlib的显示情况（默认inline是内嵌显示，通过设置为tk表示不内嵌显示）\n",
    "%matplotlib tk"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def expand(a, b, rate=0.05):\n",
    "    d = (b - a) * rate\n",
    "    return a-d, b+d"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "## 数据加载\n",
    "iris_feature = u'花萼长度', u'花萼宽度', u'花瓣长度', u'花瓣宽度'\n",
    "path = 'datas/iris.data'\n",
    "data = pd.read_csv(path, header=None)\n",
    "x_prime = data[np.arange(4)]\n",
    "y = pd.Categorical(data[4]).codes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "特征属性: 花萼长度 花萼宽度\n",
      "实际均值：\n",
      " [[ 5.006  3.418]\n",
      " [ 5.936  2.77 ]\n",
      " [ 6.588  2.974]]\n",
      "预测均值=\n",
      " [[ 5.01493896  3.4404862 ]\n",
      " [ 5.90114537  2.74385294]\n",
      " [ 6.6814044   3.0285628 ]]\n",
      "预测方差=\n",
      " [[[ 0.1194876   0.08969867]\n",
      "  [ 0.08969867  0.12147459]]\n",
      "\n",
      " [[ 0.27544608  0.08866062]\n",
      "  [ 0.08866062  0.09382524]]\n",
      "\n",
      " [[ 0.36087007  0.05158991]\n",
      "  [ 0.05158991  0.08923683]]]\n",
      "顺序：\t [0 1 2]\n",
      "准确率：79.33%\n",
      "\n",
      "\n",
      "特征属性: 花萼长度 花瓣长度\n",
      "实际均值：\n",
      " [[ 5.006  1.464]\n",
      " [ 5.936  4.26 ]\n",
      " [ 6.588  5.552]]\n",
      "预测均值=\n",
      " [[ 5.0060006   1.46399865]\n",
      " [ 6.04240777  4.41742864]\n",
      " [ 6.58888904  5.63329718]]\n",
      "预测方差=\n",
      " [[[ 0.12176525  0.01581631]\n",
      "  [ 0.01581631  0.0295045 ]]\n",
      "\n",
      " [[ 0.28119672  0.23746926]\n",
      "  [ 0.23746926  0.31503012]]\n",
      "\n",
      " [[ 0.48521779  0.36602418]\n",
      "  [ 0.36602418  0.32601109]]]\n",
      "顺序：\t [0 1 2]\n",
      "准确率：91.33%\n",
      "\n",
      "\n",
      "特征属性: 花萼长度 花瓣宽度\n",
      "实际均值：\n",
      " [[ 5.006  0.244]\n",
      " [ 5.936  1.326]\n",
      " [ 6.588  2.026]]\n",
      "预测均值=\n",
      " [[ 5.00607264  0.23754806]\n",
      " [ 6.56291563  2.02408174]\n",
      " [ 5.94928821  1.32089151]]\n",
      "预测方差=\n",
      " [[[ 0.1239802   0.01055412]\n",
      "  [ 0.01055412  0.00914172]]\n",
      "\n",
      " [[ 0.41146765  0.0558402 ]\n",
      "  [ 0.0558402   0.07689828]]\n",
      "\n",
      " [[ 0.29249903  0.07999787]\n",
      "  [ 0.07999787  0.0509237 ]]]\n",
      "顺序：\t [0 2 1]\n",
      "准确率：96.00%\n",
      "\n",
      "\n",
      "特征属性: 花萼宽度 花瓣长度\n",
      "实际均值：\n",
      " [[ 3.418  1.464]\n",
      " [ 2.77   4.26 ]\n",
      " [ 2.974  5.552]]\n",
      "预测均值=\n",
      " [[ 3.41800009  1.46400001]\n",
      " [ 2.97017899  5.56073357]\n",
      " [ 2.80062882  4.43004172]]\n",
      "预测方差=\n",
      " [[[ 0.14227691  0.01144799]\n",
      "  [ 0.01144799  0.029505  ]]\n",
      "\n",
      " [[ 0.11477629  0.07760424]\n",
      "  [ 0.07760424  0.38871245]]\n",
      "\n",
      " [[ 0.09376548  0.10702236]\n",
      "  [ 0.10702236  0.34454954]]]\n",
      "顺序：\t [0 2 1]\n",
      "准确率：92.67%\n",
      "\n",
      "\n",
      "特征属性: 花萼宽度 花瓣宽度\n",
      "实际均值：\n",
      " [[ 3.418  0.244]\n",
      " [ 2.77   1.326]\n",
      " [ 2.974  2.026]]\n",
      "预测均值=\n",
      " [[ 3.41800003  0.244     ]\n",
      " [ 2.79657869  1.31224583]\n",
      " [ 2.93629236  1.98607968]]\n",
      "预测方差=\n",
      " [[[ 0.14227697  0.011208  ]\n",
      "  [ 0.011208    0.011265  ]]\n",
      "\n",
      " [[ 0.09554395  0.04869984]\n",
      "  [ 0.04869984  0.03787478]]\n",
      "\n",
      " [[ 0.11263095  0.06192916]\n",
      "  [ 0.06192916  0.08966439]]]\n",
      "顺序：\t [0 1 2]\n",
      "准确率：93.33%\n",
      "\n",
      "\n",
      "特征属性: 花瓣长度 花瓣宽度\n",
      "实际均值：\n",
      " [[ 1.464  0.244]\n",
      " [ 4.26   1.326]\n",
      " [ 5.552  2.026]]\n",
      "预测均值=\n",
      " [[ 1.46399926  0.24399973]\n",
      " [ 4.32760641  1.36230588]\n",
      " [ 5.60395923  2.0545418 ]]\n",
      "预测方差=\n",
      " [[[ 0.02950475  0.00558391]\n",
      "  [ 0.00558391  0.01126496]]\n",
      "\n",
      " [[ 0.25410036  0.09152045]\n",
      "  [ 0.09152045  0.05088321]]\n",
      "\n",
      " [[ 0.29156635  0.03719346]\n",
      "  [ 0.03719346  0.07073705]]]\n",
      "顺序：\t [0 1 2]\n",
      "准确率：97.33%\n",
      "\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# 类别数量\n",
    "n_components = 3\n",
    "# 不同特征属性分类\n",
    "feature_pairs = [[0, 1], [0, 2], [0, 3], [1, 2], [1, 3], [2, 3]]\n",
    "plt.figure(figsize=(10, 6), facecolor='w')\n",
    "cm_light = mpl.colors.ListedColormap(['#FFA0A0', '#A0FFA0', '#A0A0FF'])\n",
    "cm_dark = mpl.colors.ListedColormap(['r', 'g', 'b'])\n",
    "\n",
    "for k, pair in enumerate(feature_pairs):\n",
    "    print (u\"特征属性:\", iris_feature[pair[0]], iris_feature[pair[1]])\n",
    "    x = x_prime[pair]\n",
    "    # 计算实际的均值\n",
    "    mean = np.array([np.mean(x[y == i], axis=0) for i in range(3)])\n",
    "    print (u\"实际均值：\\n\",mean)\n",
    "    \n",
    "    # 模型构建\n",
    "    gmm = GaussianMixture(n_components=n_components, covariance_type='full', random_state=28)\n",
    "    gmm.fit(x)\n",
    "    \n",
    "    # 效果参数输出\n",
    "    print (\"预测均值=\\n\", gmm.means_)\n",
    "    print (\"预测方差=\\n\", gmm.covariances_)\n",
    "    \n",
    "    # 获取预测结果\n",
    "    y_hat = gmm.predict(x)\n",
    "    order = pairwise_distances_argmin(mean, gmm.means_, axis=1, metric='euclidean')\n",
    "    print ('顺序：\\t', order)\n",
    "    \n",
    "    # 修改预测结果顺序\n",
    "    n_sample = y.size\n",
    "    n_types = 3\n",
    "    change = np.empty((n_types, n_sample), dtype=np.bool)\n",
    "    for i in range(n_types):\n",
    "        change[i] = y_hat == order[i]\n",
    "    for i in range(n_types):\n",
    "        y_hat[change[i]] = i\n",
    "    \n",
    "    # 计算准确率\n",
    "    acc = u'准确率：%.2f%%' % (100*np.mean(y_hat == y))\n",
    "    print (acc)\n",
    "    \n",
    "    # 画图\n",
    "    x1_min, x2_min = x.min()\n",
    "    x1_max, x2_max = x.max()\n",
    "    x1_min, x1_max = expand(x1_min, x1_max)\n",
    "    x2_min, x2_max = expand(x2_min, x2_max)\n",
    "    x1, x2 = np.mgrid[x1_min:x1_max:500j, x2_min:x2_max:500j]\n",
    "    grid_test = np.stack((x1.flat, x2.flat), axis=1)\n",
    "    grid_hat = gmm.predict(grid_test)\n",
    "\n",
    "    change = np.empty((n_types, grid_hat.size), dtype=np.bool)\n",
    "    for i in range(n_types):\n",
    "        change[i] = grid_hat == order[i]\n",
    "    for i in range(n_types):\n",
    "        grid_hat[change[i]] = i\n",
    "    grid_hat = grid_hat.reshape(x1.shape)\n",
    "    \n",
    "    # 子图\n",
    "    plt.subplot(3, 2, k+1)\n",
    "    plt.pcolormesh(x1, x2, grid_hat, cmap=cm_light)\n",
    "    plt.scatter(x[pair[0]], x[pair[1]], s=30, c=y, marker='o', cmap=cm_dark, edgecolors='k')\n",
    "    xx = 0.9 * x1_min + 0.1 * x1_max\n",
    "    yy = 0.2 * x2_min + 0.8 * x2_max\n",
    "    plt.text(xx, yy, acc, fontsize=14)\n",
    "    plt.xlim((x1_min, x1_max))\n",
    "    plt.ylim((x2_min, x2_max))\n",
    "    plt.xlabel(iris_feature[pair[0]], fontsize=14)\n",
    "    plt.ylabel(iris_feature[pair[1]], fontsize=14)\n",
    "    plt.grid()\n",
    "    print()\n",
    "    print()\n",
    "    \n",
    "plt.tight_layout(2)\n",
    "plt.suptitle(u'EM算法鸢尾花数据分类', fontsize=20)\n",
    "plt.subplots_adjust(top=0.90)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
